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. Invariant discretization schemes are derived for the one- and two-dimensional shallow-water 

f^) equations with periodic boundary conditions. While originally designed for constructing invari- 

04 ■ ant finite difference schemes, we extend the usage of difference invariants to allow constructing 

£H ■ of invariant finite volume methods as well. It is found that the classical invariant schemes con- 

verge to the Lagrangian formulation of the shallow- water equations. These schemes require to 
redistribute the grid points according to the physical fluid velocity, i.e., the mesh cannot remain 
C*") , fixed in the course of the numerical integration. Invariant Eulerian discretization schemes are 

proposed for the shallow-water equations in computational coordinates. Instead of using the 
, fluid velocity as the grid velocity, an invariant moving mesh generator is invoked in order to de- 

Qh termine the location of the grid points at the subsequent time level. The numerical conservation 

^ fin of energy, mass and momentum is evaluated for both the invariant and non-invariant schemes. 

> ■ 

a- 
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1 Introduction 

Discretization schemes that preserve characteristic properties of systems of differential equations 
QO have received increasing attention over the past years and led to the development of the field of 

geometric numerical integration. The principal motivation for this approach is that controlling 
the local discretization error, as is done in most of the classical numerical methods, can fail to 
capture essential qualitative features of the underlying problem, which might be equally impor- 
■ tant in order to obtain reasonable integration results. Such features can include, but are not 

CN ■ necessarily limited to, conservation laws, point symmetries, Hamiltonian structure, conservation 

of phase-space volume and asymptotic characteristics. Various geometric numerical integration 
schemes have been developed that capture these properties in the course of discretization, such as 
^ ■ conservation laws and the Hamiltonian structure [10, 34, 43], Lie symmetries [3, 20, 21, 22, 33, 50] 

^ . and phase-space volume [26, 47]. 

In the present paper, we aim to concentrate on the problem of deriving discretization schemes 
with symmetry properties by developing appropriate finite difference and finite volume schemes 
for the shallow- water equations. In particular, we are concerned with the problem of finding 
discretization schemes that are invariant under the maximal Lie invariance group admitted by 
the shallow-water equations with (double) periodic boundary conditions. Choosing the shallow- 
water equations for such an investigation can be motivated because they constitute a prominent, 
simple, yet fully-nonlinear model of fluid mechanics exhibiting various features of the original 
set of governing equations of hydrodynamics, such as the simultaneous occurrence of both fast 
(divergent) and slow (vortical) waves and the existence of conservation laws, symmetries and 
a Hamiltonian form. Moreover, the shallow-water equations always served as an important 
intermediate model to test new numerical schemes [2, 28, 43, 44, 45, 46]. 

Similarly to conservation laws, symmetries have important implications on the solutions of 
differential equations. When simulating the dynamics of a classical mechanical system in a con- 
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stantly moving coordinate system, it should be a clear desire that the numerical model to be 
used for that problem is Galilean invariant as otherwise physical laws can be violated. It is also 
well known [15, 16, 29] that the shape of a solution of a system of differential equations near 
a blow-up point can tend to a group-invariant solution of this system. Often group-invariant 
solutions well capture so-called intermediate asymptotic behavior of the solutions after a suffi- 
ciently long period of evolution. For the simulation of invariant solutions, symmetry-preserving 
discretization schemes can give better numerical results than standard schemes that do not 
preserve the geometry of differential equations. 

The design of invariant discretization schemes for evolution equations in general requires the 
explicit treatment of meshes that are not time-space orthogonal, i.e., time-adaptive grids. Such 
grids pose several challenges from the numerical point of view that have not been well investi- 
gated in the field of invariant numerical schemes up to now. On the other hand, meshes that 
adapt according to the development of the numerical solution are an extensively investigated 
subject in the field of numerical mathematics, see, e.g., [29, 49]. The question not explicitly 
answered so far is whether the problem of finding discretization schemes with symmetry prop- 
erties can be embedded into the study of adaptive numerical schemes in the multidimensional 
case. In the present paper, we discuss a possible answer to that problem, exemplified by the 
shallow- water equations. 

The outline of the paper is as follows: Properties of the shallow-water equations are discussed 
in Section 2. Section 3 is devoted to a review of common techniques that allow one to construct 
invariant finite difference schemes. In Section 4 we derive invariant discretization schemes for the 
one-dimensional shallow-water equations. This is done both by using the Lagrangian description 
of the shallow-water equations and by setting up an invariant grid generator for Eulerian schemes 
in computational coordinates. In Section 5 we discuss strategies for the design of invariant 
numerical models in higher dimensions and illustrate them with the two-dimensional shallow- 
water equations. Again, both Lagrangian schemes in physical coordinates and Eulerian schemes 
in computational coordinates with an invariant grid generator are introduced. For the first 
scheme we use an invariant finite volume discretization, while the second scheme is based on 
finite differences. A summary and concluding remarks can be found in Section 6. 



2 Symmetries of the shallow-water equations 

The nondimensionalized system of shallow-water equations in Cartesian coordinates is 
u t + uu x + vuy + h x = 0, 

v t + uv x + VVy + h y = 0, (1) 

h t + uh x + vh y + h(u x + v y ) = 0, 

where v = (u, v) is the fluid velocity in the plane and h is the height of the fluid column over 
a fixed reference level within the fluid. The bottom topography is assumed to be flat here for 
simplicity. Treating non-flat topographies would lead to the inclusion of additional source terms 
in system (1). The shallow- water equations are derived from the Euler equations for an ideal 
fluid under the following assumptions: the validity of the hydrostatic approximation, constancy 
of the fluid density and much smaller scale of vertical motions in comparison with horizontal 
motions [41]. 

The shallow- water equations (1) can be represented in Hamiltonian form [35] using 

as a Poisson bracket, where J- and Q are functionals of v and h, q = Q/h = (v x — u y )/h is the 
potential vorticity, k denotes the vertical unit vector, &A = dxdy is the area element, and the 
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integration extends over the domain of the entire fluid. The Hamiltonian for the shallow-water 
equations is given by the total energy 



'H = \J {hv 2 + h 2 ) dA. 



Additional conserved quantities are associated with the above noncanonical Poisson bracket. 
For any function / of the potential vorticity q, the integral 



Cf = j hf(q)dA 



is conserved on solutions of the shallow- water equations. This class of conserved quantities 
contains the mass M = C\ , the circulation Z = C q and the potential enstrophy 6 = C q 2 / 2 . Two 
more conserved quantities are the momenta in the x- and y-directions, 



V x = J hu dA, V y = J hv dA. 

The maximal Lie invariance algebra 02 of the two-dimensional shallow- water equations (1) is 
nine-dimensional; see, e.g., [19, 40]. A basis of this algebra consists of the vector fields 

d t , d x , d y , td x + d u , td y + d v , 

td t + xd x + yd y , xd x + yd y + ud u + vd v + 2hd h , 

— yd x + xd y — vd u + ud v , t 2 dt + txd x + tyd y + (x — tu)d u + (y — tv)d v — 2thdh- 

These vector fields generate one-parameter Lie symmetry groups, which correspond to (i) time 
translations, (ii)-(iii) space translations, (iv)-(v) Galilean transformations, (vi)-(vii) scalings, 
(viii) rotations and (ix) inversions in t. 

In what follows, we will also use the one-dimensional version of system (1), in which case we 
set v = and drop the dependence of u and h on y. The resulting system reads 

u t + uu x + h x = 0, h t + uh x + hu x = (2) 

and preserves the one-dimensional versions of total energy, mass and momentum, 



n = ^ J (hu 2 + h 2 ) dx, M = j hdx, V = j hudx. 



It is well known that the maximal Lie invariance algebra gi of system (2) is infinite dimen- 
sional and spanned by the vector fields 

td t + xd x , xd x + ud u + 2hd h , td x + d u , 

(2s - 6tu)d t + (6/1 - 3u 2 )td x + [u 2 + 4h)d u + 4hud h , f{h, u)d t + g{h, u)d x , 

where the functions / and g run through the set of solutions of the system 

9h ~ uf h + f u = 0, g u ~ uf u + hf h = 0. (3) 

The existence of the latter generator is owed to the possibility of linearization of system (2) to 
system (3) by means of the hodograph transformation in that u and h are assumed as the new 
independent variables and / = t and g = x are the new unknown functions. The linearization 
by the hodograph transformation permuting the pairs of dependent and independent variables 
is a general property of homogeneous first-order systems of partial differential equations in two 
independent variables and two unknown functions that are linear in derivatives with coefficients 
depending only on the unknown functions. See also [30, p. 154] for the symmetry interpretation of 
the linearization of the one-dimensional shallow-water equations. Note that system (3) is reduced 
to a single Tricomi equation. More precisely, excluding g by cross differentiation, we obtain the 
equation f uu = hfhh + 2/^. The substitution <f> = hf then leads to the Tricomi equation 
4>uu = hcfthh- Another way to reduce the system is to rewrite the equation f uu = hfhh + 2fh m 
the form h 3 f uu = h 2 (h 2 fh)h an d to carry out the transformation z = 1/h, which yields a similar 
Tricomi equation, f uu = z 3 f zz . Symmetry analysis of such equations was carried out in [7, 9]. 
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3 Construction of invariant numerical discretization schemes 



The problem of constructing discretization schemes that preserve symmetries of the correspond- 
ing differential equations was first systematically addressed by Dorodnitsyn and his collabora- 
tors [3, 11, 22, 23, 24, 25]. As there are an infinite number of possibilities to approximate a differ- 
ential equation by means of finite differences, one might single out those among all possible dif- 
ference schemes that inherit symmetries of the original differential equation. Dorodnitsyn's ap- 
proach can be summarized in the following way: First determine the maximal Lie invariance alge- 
bra of the model under consideration. For many classical hydrodynamical problems, this task has 
already been completed [31, 39]. Then a discretization stencil has to be chosen. The generators of 
Lie symmetries are then prolonged to all points of the stencil. From these prolonged generators, 
the invariants of the extended group action on the stencil are determined. The final step is then to 
assemble the obtained invariants into a difference approximation of the original differential equa- 
tion. By difference approximation it is meant that in the continuous limit, the invariant finite dif- 
ference scheme reduces to the original differential equation in some coordinates. Each difference 
approximation consists of a physical difference equation and equations governing the positions of 
grid points. In the continuous limits, these grid equations often reduce to some trivial identities. 

Altogether, this method is a straightforward application of inverse group classification, using 
transformation groups acting on functions defined on a discrete set of points rather than on 
a continuous space. In the usual inverse group classification one starts with a particular Lie 
group G and aims at finding those systems of differential equations that admit G as a symmetry 
group. In practice these systems are found by computing differential invariants (i.e., invariants 
that involve derivatives of dependent variables) of G. Any function of differential invariants is a 
differential invariant of G and, subject to some regularity condition, any system of differential 
equations can be expressed in terms of differential invariants of its maximal Lie invariance 
group [38]. The Dorodnitsyn method works by selecting the maximal Lie symmetry group of a 
system of differential equations as the initial Lie group G. By extending the action of G to the 
points of the discretization stencil, one is able to compute invariants of the extended action, i.e., 
difference invariants of G. As in the continuous case, any function of difference invariants is a 
difference invariant. Constructing a difference approximation of a system of differential equations 
using difference invariants therefore leads to a symmetry-preserving discretization scheme. 

A common feature of difference schemes constructed by the above method is that grid points 
might not remain fixed in the course of the numerical integration. Precise criteria for a grid 
to be uniform, orthogonal and possessing flat time layers are formulated as conditions on the 
coefficients of infinitesimal symmetry generators and are broken, e.g., for Galilean boosts and 
inversions [22]. This means that for such symmetries it is not possible to use isotropic or static 
grids. Hence, the problem of establishing good conditions governing the position of grid points 
both spatially and temporally becomes vital. 

Up to now the reviewed technique has been applied to physically rather simple models, usually 
only involving time and one space dimension [3, 11, 24, 25, 50]. It is understandable that the 
multidimensional case is even more delicate, as there is an increasing number of possibilities for 
assembling the difference invariants to finite difference schemes. In addition, grids can evolve 
differently in distinct spatial dimensions, which might cause severe numerical problems, such as 
tangling meshes, if not treated appropriately. In the present paper, we aim to discuss ways of 
overcoming the latter problem. 

An alternative approach to constructing finite difference schemes with symmetry properties 
uses moving frames in the Fels and Olver formulation [27]. In contrast to the Dorodnitsyn 
approach, where finite difference schemes are constructed from the outset, in the moving frame 
method the concept of invariantization of existing schemes plays the key role. This technique 
can be summed up as follows [32, 33]: Determine the Lie symmetry group of the given system of 
differential equations. This part is standard and usually involves exponentiation of elements of 
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the maximal Lie invariance algebra of the system. Subsequently, a moving frame associated with 
the Lie symmetry group is constructed. Roughly speaking, a moving frame is an equivariant 
function that returns the unique group element mapping a given point to a point of a chosen 
submanifold (the cross-section), which intersects each group orbit once and transversally. Since 
the condition for submanifolds to be cross-sections to the group orbits is quite general, there 
is a freedom in choosing such submanifolds and hence in constructing the associated moving 
frames. Once a moving frame is obtained, it can be used to map an arbitrary function to an 
invariant function. This is a general property of any moving frame that allows determining 
usual invariants and differential invariants of a group action [36, 37]. In the same way, it is 
possible to take a given difference scheme (considered as a function of grid points) for a system 
of differential equations and invariantize it by a moving frame. By this procedure, the given 
scheme is transformed to a new scheme that will be invariant under the same Lie group as was 
used to determine the moving frame. 

The main benefit of this method is the possibility of using existing finite difference schemes as 
a starting point in the development of invariant schemes. Consequently, such invariant schemes 
could eventually be implemented in existing numerical models with limited effort. At the same 
time, the freedom in constructing a moving frame can make it somewhat difficult to predict the 
precise form of invariantized expressions and, therefore, to arrive at a scheme that not only is 
invariant but also has some desirable numerical properties, as discussed, e.g., in [20]. Although 
assembling difference invariants to invariant discretization schemes can be realized in a variety of 
ways too, the form of the particular difference invariants usually imposes enough hints in order 
to find reasonable finite difference approximations of a given system of differential equations. As 
with the Dorodnitsyn method described above, invariantization of existing numerical schemes 
may also lead to grids that evolve during numerical integration. 

Another method for the construction of schemes with certain invariance properties was pro- 
posed in [14] for equations describing blow-up problems; see also [15, 16, 29]. The main idea in 
this approach is to use adaptive moving meshes from the very beginning because they are well 
suited for problems that develop shocks after a finite integration time. As a moving mesh compli- 
cates the discretization of differential equations in the physical space, the system to be discretized 
is first transformed into so-called computational coordinates that remain orthogonal and do not 
evolve during the numerical integration. The physical system is then discretized in the compu- 
tational coordinates. It is advocated in this approach that for equations exhibiting blow-ups and 
for the description of the solution near the singularity, scale invariance plays an exceptional role. 
Therefore, scale invariance is required to be preserved in the course of discretization. The evolu- 
tion of the mesh is formulated as an auxiliary system of differential equations, the so-called mov- 
ing mesh partial differential equations. The auxiliary system is then selected in such a manner as 
to preserve the scale invariance of the original physical model. A straightforward extension to the 
above approach is to require the mesh equations to possess not only the scale symmetries but also 
the other symmetries that the system of physical differential equations admits. The following 
property is basic for this extension: The prolongation of any point symmetry of the initial sys- 
tem L of differential equations to the computational coordinates by means of the identical trans- 
formation is a point symmetry of the counterpart of C in terms of the computational coordinates. 

In the present paper, we will introduce yet another approach to the construction of invari- 
ant discretization schemes, which will be essential for multidimensional systems of differential 
equations. It rests on first expressing the system of differential equations under consideration 
in terms of computational coordinates and then extending the symmetry transformations of the 
original system to the system written in computational variables. Once it is understood how 
the system behaves under the extended symmetry transformations, one constructs a finite dif- 
ference scheme that is transformed by the discretized version of the extended transformations in 
a similar way. In addition, the extra differential equations that control the location of the grid 
points are discretized in an invariant way, e.g., by using the finite difference invariants. 
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The reason why it is necessary to develop one more technique for the construction of invariant 
discretization schemes is twofold. Firstly, it is rather difficult to set up a proper invariant 
scheme for systems of differential equations using difference invariants as basic building blocks. 
For multidimensional moving meshes, an additional problem is to find proper finite difference 
analogs of derivatives. Secondly, it can be (and, in general, will be) desirable to include additional 
qualitative properties of differential equations in the construction of invariant discretization 
schemes. Within the invariantization technique it might be tedious to ensure the numerical 
preservation of certain conservation laws, even if the initial system includes equations represented 
as conservation laws, which is precisely the case for the shallow-water equations written in 
momentum form. A similar remark holds for the Dorodnitsyn approach. An exception is given 
for equations derived from a variational principle for which, in view of the discrete version of 
the Noether identity, conservation laws and symmetries can be simultaneously preserved in the 
course of a proper invariant discretization of the associated Lagrangian [11]. There is still no 
algorithm using only difference invariants (or the invariantization map) that guarantees that the 
resulting invariant scheme will admit certain conservation laws. On the other hand, with the 
new approach to be introduced in the present paper, it is possible to construct, in a quite direct 
way, schemes that both are invariant and preserve some of the conservation laws possessed by 
the initial system. 

4 Invariant numerical models for the one-dimensional 
shallow-water equations 

In Section 2 we discussed the Lie symmetries of the shallow- water equations without any relation 
to boundary value problems. However, when setting up a numerical model for a specific set of 
problems, the explicit treatment of certain boundary conditions is usually inevitable. As a rule, 
the Lie symmetries possessed by a boundary value problem form only a subgroup (often even 
trivial one) of the maximal Lie symmetry group admitted by the involved system of differential 
equations [8]. Stated in another way, a specific boundary value problem usually admits only 
a small subset of the symmetries of the associated system of differential equations considered 
without boundary and initial conditions. This is in particular the case for various differential 
equations arising in hydrodynamics, which admit wide Lie invariance groups in the absence of 
boundary and initial conditions [1, 4, 5, 6, 31, 39]. As shown in the present section, this is also 
the case for the shallow-water equations. Therefore, it is necessary to find a way to incorporate 
the boundary and initial conditions considered into the numerical model to be developed. 

4.1 Selection of symmetries using boundary conditions 

In order to design invariant numerical schemes for a system of differential equations, two principal 
strategies can be adopted. 

In the first approach, which is applied in most of the previous works on invariant discretiza- 
tion [22], numerical schemes preserving the entire maximal Lie invariance groups of the corre- 
sponding systems are developed and then implemented for the specific physical configurations 
of interest. The drawback of this approach is that the practical implementation of a numerical 
scheme always requires the explicit treatment of a boundary value problem. As was said above, 
for the boundary conditions arising most often in hydrodynamics (e.g., periodic, reflective or ab- 
sorbing) , the maximal Lie invariance groups of systems without boundary conditions are usually 
much wider than those of particular boundary value problems. Using the first approach may 
therefore lead to the overly restrictive requirement that all the symmetries of the considered 
system of differential equations are equally important in the course of invariant discretization. 
For a particular boundary value problem, however, this may not be the case, as some of the 
symmetries of the system might not be admitted at all. 
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This is why we adopt the second approach here, which only requires the preservation of 
symmetries that are compatible with the class of specific boundary value problems under con- 
sideration. The apparent drawback of this approach is that, if one aims to test different kinds 
of boundary conditions, it can be necessary to design a new scheme for each configuration, as 
different symmetry groups may arise when varying specific settings. On the other hand, for a 
model to be used for a particular purpose (e.g., a weather or climate prediction model), the 
boundary conditions are generally fixed at the stage of model development and therefore do not 
change subsequently. Another advantage of the second approach is of more physical nature. The 
restriction imposed on symmetries in that they map a given boundary value problem to itself 
is unreasonable even from the physical perspective. When transforming a given reference frame 
to another reference frame, the boundary value problems of the reference frames involved are 
also mapped to each other. Therefore, it is not natural to require appropriate symmetries of a 
system of differential equations to preserve a particular boundary value problem but rather only 
to impose that these symmetries map boundary value problems from a class of such problems 
(e.g., periodic domains of any size with varying initial time and initial conditions) to each other. 
Such transformations are known as equivalence transformations, and when deriving symmetry- 
preserving discretization schemes, we require a subgroup of the maximal Lie invariance group 
of the original system of differential equations to be compatible with the structure of a prede- 
fined class of boundary value problems, i.e., elements of the subgroup should act as equivalence 
transformations on the boundary conditions rather than as symmetry transformations. 

The relaxed condition of requiring the finite difference schemes to be invariant only under the 
transformations admitted by a class of boundary value problems thus provides a natural selection 
criterion for subgroups of a maximal Lie invariance group to be preserved numerically. In view 
of the particular nature of the infinite-dimensional maximal Lie invariance algebra 51 of the 
one-dimensional shallow-water equations (or, more generally, systems of differential equations 
arising in hydrodynamics), it might be a cumbersome or even useless task to attempt to preserve 
all these symmetries in the respective discrete models. 

Among the most natural boundary conditions for the one-dimensional shallow-water equa- 
tions in the setting of geophysical fluid dynamics are periodic ones, which we aim to study here. 
These boundary conditions are advantageous as they are not as restrictive as, e.g., Dirichlet 
boundary conditions from the pure symmetry point of view and generally lead to the selec- 
tion of symmetries that have a clear physical interpretation. The subalgebra Si of Qi that is 
compatible with periodic boundary conditions is spanned by the vector fields 

d t , d x , td x + d u , td t +xd x , xd x + ud u + 2hd h . (4) 

Even if we neglected initial conditions, only elements from the narrower subalgebra (dt,d x , 
td x + d u , tdt — ud u — 2hdf l ) generate one-parameter symmetry groups of such a boundary value 
problem, as scalings with respect to x are not admitted once a domain (periodicity) length is 
fixed. However, the inclusion of scaling symmetries in the subsequent consideration is justified 
as they are equivalence transformations of the chosen class of boundary value problems. In other 
words, upon preserving the subalgebra Si we are still able to test different domain lengths. In 
this sense, the preservation of scaling transformations plays an important role for the class of 
boundary value problems for the one-dimensional shallow-water equations with periodic bound- 
ary conditions and any domain size even if scalings are not proper symmetry transformations 
for a specific numerical integration. 

4.2 Classical invariant schemes and beyond 

We begin our study of invariant numerical schemes for the shallow-water equations using the 
classical construction proposed by Dorodnitsyn. Within this framework, we have to prolong 
the selected subalgebra Si to the discretization stencils that we aim to use. These stencils are 
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depicted in Fig. 1. As the symmetry group associated with Si does not violate the criterion for 
using flat time layers (see [22] for more details), all points in the spatial domain are defined at 
the same time. However, if we wish to preserve Galilean invariance in a numerical scheme, it is 
impossible to use a fixed grid, i.e., X{ ^ Xj. Here and in what follows variables with a hat and 
without a hat denote values on the grid at the time levels t + r and t, respectively, and r is the 
time step. The possibility of evolving grids in general also leads to nonhomogeneous spacings in 
the course of the integration; i.e., xi + \ — xi ^ &i — even if the initial grid {x^ is equally 
spaced. 



(i + r, ii-i , Uj_i, hi-i) (t + T,Xi,Ui,hi) 



(i + r, x i+1 ,u l+1 ,h l+1 ) 



(t,Xi,Ui,hi) 



(t,x l+1 ,u, +1 ,h r+ i) 



Figure 1. Stencils for invariant schemes of the one-dimensional shallow-water equations. An explicit (Euler 
forward) and an implicit (Euler backward) scheme is defined using the points indicated by filled and dashed 
circles, respectively. 

Prolonging the vector fields (4) to the points indicated in Fig. 1 gives 

d t , d Xi + d Xi+1 + 8 Xi _ x + d Xi + d Xi+1 + , 

t(d x . + d Xi+1 + d Xi _,) + (t + r){d Xi + % +1 + ^._J + 
d Ul + d Ui+1 + cV, + du, + du i+1 + , 

td t + rd T + Xid Xx + x i+1 d Xi+1 + Xi-\d Xi _^ + Xid x . + x i+1 d Xi+1 + Xi-\d Xi _ x , (5) 

Xid x . + x i+1 d Xi+1 + + Xid Xi + x i+1 d x . +1 + x^d^ + 

uid Ui + u i+1 d Ui+1 + Ui-A^ + Uidu r + Ui+idn i+1 + Ui-xdn^ + 

+ 2h i+1 d hi+1 + 2^_i^ i _ 1 + 2/t^. + 2h i+ id h . +l + 2^-1^. 

To construct an explicit (Euler forward) numerical scheme, we have to restrict ourselves in (5) 
to the values that are defined at the filled circles depicted in Fig. 1. A convenient complete set 
of functionally independent difference invariants is then 



lo = — , h 

U i+ i - Ui-i 

3 = — z — 

hi 

(X i+ l - 



Xi Hi 



-r, h 



r, I 



r, /. 5 



Uj-Uj 

3<i— 1 



2 T 



(6) 



(Zj+l - Xj_l) 5 



(X i+ l - Xi-lf 



where Xi = (xi — Xi)/r is by definition the mesh velocity. These invariants are found from 
integrating the system of first-order quasilinear partial differential equations Vj(J) = 0, where 
Vj, j = 1, ... ,5, are the prolonged vector fields presented in (5). The determining system for 
invariants admits precisely m\ = m s — r functionally independent solutions, where m s is the 
number of stencil variables and r is the rank of involved vector fields. For the explicit scheme, 
we have m s = 14, r = 5 and hence m; = 9. 
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Using the difference invariants of the set (6), we can approximate (2) via 



h = o, i 2 + i 7 -h = o, h-h + hh = o 

or, explicitly, 

Ui-Ui h i+1 - hi-i hi - hi u i+ i - Ui-i 
Xi = Ui, 1 = 0, \-hi = 0. {<) 

T X i+ i - Xi-l T X i+ l - Xi-l 

In the continuous limit the above scheme leads to the following system of differential equations: 
dx du dh dh , du 

d! = U < dt + &=°' d? + /i ^ = °' (8) 

which is (2) in Lagrangian variables. 

As the Euler forward scheme is only conditionally stable, it is beneficial to construct an 
implicit invariant numerical scheme. A simple implicit scheme is the Euler backward scheme, 
which can be constructed in a similar way as the invariant Euler forward scheme. However, we 
prefer to at once construct a trapezoidal scheme, which has in general a greater accuracy. To 
accomplish this we additionally need the difference invariants 

±i — Ui Ui + i — hi + i — 2 , _» 

-t, ho = z r, In = - — ; -t\ (9) 



-Xi-l X i+ l-Xi-l (X i+1 - Xi-l)[X i+ l - Xi-l) 

Note that {Iq, . . . , In} is not a complete set of functionally independent invariants for the 
transformation group generated by the vector fields (5) on the trapezoidal stencil as the total 
number of variables on this stencil equals 20 and hence a functional basis of related invariants 
consists of 15 invariants. By combining the invariants and /5-/11 we construct 

h+I 9 = 0, h + \{h-h + hi) = Q, I 8 -h + ^(hh + hho)=0, 
which boils down to the form 

4-1 — hi— 1 /l,'4-1 — \ 

(10) 

0. 

Xi+l Xi—% / 

This scheme also converges to (8). 

The problem with the above schemes in particular and with the difference invariants ap- 
proach to the construction of invariant schemes in general is that it is hard to control properties 
other than symmetries that the resulting discretizations admit. It can be checked by direct 
computation that the above schemes violate even the mass conservation law (conservation of 
momentum and energy is violated as well). This violation of fundamental conservation laws 
is a direct consequence of the construction method of the invariant finite difference schemes, 
which only takes into account local information on itj and hi (i.e., the difference invariants) but 
provides no guideline ensuring the preservation of global features by the numerical solution. 

In the present case, this problem can be partially circumvented by discretizing the shallow- 
water equations not in the form (8) but rather in the momentum form 

Eq h =h t + {uh) x = 0, Eq u = (uh) t + (hu 2 + h 2 ^j = 0. (11) 
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An invariant finite difference approximation of this system using an Euler time step is 



i Xi+i - Xi-x ~ x i+x - Xi-i t h 2 i+l - h\_ x 
Xi = Ui, hi hi = 0, Uihi Uihi + - = 0, 

while the invariant trapezoidal scheme is 




(12) 



Uihi mhi + - — h — 11 = 0. 



The schemes constructed in this way numerically preserve the mass conservation law. The 
explicit scheme in addition preserves the momentum, while the implicit scheme (12) does not 
preserve momentum exactly, as — ^ Xi + \ — x%—\ in general. This condition, however, 
would be required to yield exact conservation of momentum but, as the change in the spacings 
Xi + \ — Xi-i is not abrupt, the violation of momentum conservation of the scheme (12) is rather 
small; see the result of a numerical integration using (12) below. None of these schemes respects 
conservation of energy. 

In the continuous limits both the explicit and implicit discretizations give 

dx dh du d(uh) 1 dh 2 

d^ = n ' Tt +h d- x =^ ^T + 2^ = ' 

which is (11) expressed in Lagrangian variables. Similar as the schemes (7) and (10) the above 
two invariant schemes for the momentum form of the shallow-water equations could be expressed 
in terms of difference invariants. As these expressions are considerably more involved than the 
analogous expressions for (7) and (10), we do not present these difference invariant forms here. 
The reason for the difference invariant expressions being more complicated in the present case 
can be traced back to the result of the Galilean transformation t = t, x = x + e±t, u = u + e±, 
when applied to (11), which yields 

Ec> = Eq h , Eq u = Eq u + £iEq h . (13) 

The transformation of the momentum equation thus yields a linear combination of the momen- 
tum equation with the continuity equation. This implies that the momentum equation itself 
cannot be expressed in terms of differential invariants but only in combination with the conti- 
nuity equation. It is thus not natural to approximate the momentum equation using difference 
invariants. At the same time, checking that the proposed conservative schemes are indeed invari- 
ant can be shown directly by acting with the prolonged vector fields (5) on them and verifying 
that the results of these operations yield zero on the solutions of the numerical scheme. 

The result of a numerical integration taking harmonic initial conditions for both u and h 
using the scheme (12) is depicted in Fig. 2. As the evolution of the mesh points is directly 
coupled to the (initially harmonic) physical velocity, the single mesh points quasi-oscillate around 
their initial positions (Fig. 2d). No special ability of the mesh to follow the developing shock 
(Fig. 26, showing the numerical solution of h at time t = 3) is visible, which is one of the major 
disadvantages of the scheme (12). The scheme conserves mass up to machine precision (10~ 16 ) 
but it dissipates energy, with the relative change (H(t) — H(0)) fH(0) of energy being of the 
order 10 -5 at the end of the integration. The relative change in momentum in this integration 
is of order 10~ 14 without a positive or negative trend. The values of A4, H and V at time t 
are evaluated using the formulas M. = \ J2i hi( x i+i — Xi-i), Ti = \ ^i{hiuf + hf){xi + \ — Xj_i) 
and V=2 ^2i hiUi{xi + \ — Xj_i), respectively. 
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Figure 2. Numerical integration of the one-dimensional shallow-water equations (2) using the scheme (10) with 
r = 0.001 and N = 51 grid points on the domain [0, 2n] over the time interval [0,3]. The initial conditions are 
u = Asmx and h = ho + Asin(a; + <po), with A = 0.4, ipo = n/6 and ho = 10. (a) Evolution of the discretization 
grid. (b) Numerical solution for h at t = 3. 

4.3 Invariant discretization on equidistributing meshes 

So far, we have mainly been concerned with assembling difference invariants in a proper way, 
so as to guarantee the invariance of the resulting finite difference schemes. That is, the in- 
variance condition was the relevant starting point in the design of the above schemes. The 
main problem with this approach is the lack of an explicit error control for the proposed nu- 
merical models. When setting up a numerical scheme, as well as being of primary interest to 
ensure the discrete preservation of qualitative features of differential equations, it is also im- 
portant to address classical issues from numerical analysis. For adaptive moving meshes, these 
issues mainly concern the prevention of abruptly changing grids, mesh racing and mesh tan- 
gling, which can significantly degrade the numerical solution and ultimately lead to convergence 
failure. 

When dealing with finite difference schemes on adaptive moving meshes, one usually re- 
gards the mesh movement as a time-dependent coordinate transformation from a fixed logical 
(computational) domain to the physical domain of the system of differential equations. The 
computational coordinates are defined to index the positions of the grid points in the mesh. 
Because in any regular grid each grid point keeps its position relative to its neighbors in the 
mesh even in the presence of adaptation, it is convenient to take the (spatial) computational 
coordinates as time-independent, Cartesian and orthogonal, with uniform spacing on the unit 
interval (up to scaling). In the one-dimensional case considered here, the step of the spatial 
computational coordinate £ equals 1/(N — 1), where N is the number of grid points at each 
fixed time level. In order to use computational coordinates, it is necessary to transform the 
system of differential equations from the physical space to the index space; see, e.g., [15, 29]. 

The relation of the usage of computational coordinates to invariant numerical schemes will 
be illustrated again with the one-dimensional system of shallow- water equations. The central 
idea is that any finite difference discretization of the shallow-water equations on a moving mesh 
in computational coordinates is invariant under the Lie group generated by the vector fields (5). 
Indeed, under the transformation t = 9, x = x(9,£) to the computational coordinates (9,£), the 
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one-dimensional shallow-water equations (2) take the form 

ue + (u — xq)^- + — he = 0, h e + (u - x e ) — + — hut = 0, (14) 

X£ X£ X£ X£ 

where u = u(0,x(9,£)) and h = h(9,x(0,£)). It is obvious that any usual finite difference 
discretization of (14) possesses the symmetry group requested. For example, discretizing using 
forward differences in time and central differences in space, from the above system we obtain 

Ui-Ui u i+ i - Ui-i hi+i - hi-i 
h [Ui - Xi) 1 = 0, 

T X i+ l - Xi-l X i+ i - Xi-l 

hi — hi hi + \ — hi-x Ui + \ — Ui-i 
h (Ui - Xi) h hi = 0, 

T X i+ l - X i+ l - Xi^l 

where Uj = u(0,£i) = u(t,Xi(t)), hi = h(6,£i) = h(t,Xi(t)) and in and hi denote the same values 
at 9 + t. This discretization coincides with the second and third equations of the system (7) if 
we assume the grid evolution to be Lagrangian of the form = itj. 

In much the same way, an invariant implicit discretization (trapezoidal rule) can be obtained 
from (14), giving 

'u i+ l - Ui-l Uj+l - Ui-l \ 
Xi+l ~ Xi-l J 




0, 



Xi+l Xi—i Xi+l Xi—i 



hi + l f Ui + Uj _\ I hj+i - hj-x + hj+i - hj-i \ + 



t 2 V 2 J \ Xi+i - Xi-i Xi+i - Xi-i 



2 \ Xi+i -Xi-i Xi+i -Xi-i) 

Again, in the Lagrangian case Xi = (uj + Ui)/2, this scheme coincides with the scheme (10). 

Neither of these two schemes preserves mass and momentum, as they approximate the rep- 
resentation (14) of the shallow-water equations, where the equations are not in conserved form. 
It is possible to discretize the conserved form (11) using computational variables as well, which 
boils down to 

~ h Xi+l -Xi-l _ ^ _ rA , h , = ^ ^ Xi+l - Xi-l _ ^ + + T 2 = Q 

Xi+i - Xi-i Xi+i - Xi-i 2 

for the explicit Euler scheme (preserving mass and momentum) and to 



(15) 



k Xl+1 Xi ~ l -h- T -{A(h) + A{h)) = 0, 
Xi+i ~ Xi-i 2 

Ui hi Xl+l ~ Xi ~ l ~ Uihi + T -(A{uh) + A(uh)) + ^-(D(h 2 ) + t)(h 2 )) = 
Xi+i - Xi-i 2 4 

for the implicit trapezoidal discretization. In both schemes we denote 

At \ - ( Ui+l ~ ~ ~ x i-i) z i-i nt \ - Zi+l ~ Zi ~ l 
A\z) — — , iJ(z) — — , 

Xi+l Xi — l Xi+l Xi—l 

and A(z) and D(z) have the same forms as A{z) and D(z) with all the variables replaced by 
the associated variables on the next time step 9 + t, only keeping the grid velocity the same. In 
the continuous limit, these schemes converge to 

(x c F*) + (F x - F'xe)^ = 0, 
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which is indeed (11) in computational coordinates using F* = (h, hu) and F x = (hu, hu 2 + \h 2 ). 

The above observation can be easily extended to other invariant schemes for evolution equa- 
tions admitting Galilean transformations as symmetries. Its main benefit is that it allows us 
to establish a connection to the theory of discretization on adaptive moving meshes. This may 
aid in tackling the problem of finding invariant finite difference schemes that also have good 
numerical properties. 

In order to complete the invariant schemes in computational coordinates, it is necessary to 
determine the mesh velocity Xi in an invariant way. This can be done using equidistributing 
meshes. Classically, a mesh is called equidistributed if the relation 

x(0 rb 

p(x)dx = £ / p(x)dx 

J a 

holds for the continuous mapping x = : [0,1] — > [a,b]; see, e.g., [29]. The function p = 
p(x) is called the monitor function. It determines the regions of concentration of the grid. 
Differentiating this relation twice with respect to £, one obtains the equation 

{pxfa = 0, (16) 

with the boundary conditions x(0) = a and x(l) = b, which is satisfied for an equidistributed 
mesh. As we are considering periodic boundary conditions, we should modify the classical 
framework of equidistributing meshes and replace the boundary conditions for x(£) by setting 
x(l) — x(0) = 2tt and x^(0) = x^(l). The periodic conditions for agree with the invariance 
requested. 

The above schemes in computational coordinates will therefore be completely invariant if we 
obtain the grid on the next level (and therefore the grid velocity Xi = (xi — Xi)/r) from an 
invariant discretization of the equidistribution principle (16). The discretization 

(pi+i + Pi)(x i+i - Xi) - (pi + Pi-i){xi - = (17) 

is invariant provided that we choose an invariant monitor function p. An ansatz for p motivated 
from the theory of adaptive grids is, e.g., the arc- length (-like) monitor function 



au 



2 




with a being the (positive) adaptation constant. This monitor function is invariant with respect 
to vector fields (4) excluding only the scale operator tdt + xd x but the corresponding scalings 
are equivalence transformations for the set of such monitor functions, where the parameter a 
varies. The above ansatz for p can be discretized in an invariant way via 

(18) 

The resulting form of (17) can then be solved either using an iterative method, such as Jacobi 
or Gaufi-Seidel iteration, or by relaxation, e.g., using the moving mesh PDE approach [15, 29]. 

Remark 1. For the equation (16) to possess a Lie symmetry algebra q that is contained in the 
linear span s\ of vector fields (4) trivially extended to £, it suffices for the monitor function p to 
be an invariant of g. On solutions of the shallow-water equations (2) we can assume without loss 
of generality that the function p does not depend on derivatives of u and h involving differentia- 
tion with respect to t. Then the general form of p that is an invariant of the pure Galilean algebra 
(dt, d x ,td x + d u ) is given by an arbitrary smooth function of derivatives of u and h with respect 
to x including h itself but not u. In order to attain invariance with respect to scale transforma- 
tions, the function p should depend only on specific products of powers of the above derivatives. 
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At the same time, the incorporation of geometric properties of solutions (e.g., the length of a 
graph between neighboring grid points) into the monitor function is more important than scale 
invariance. Therefore, scale transformations can be allowed to act in a relaxed way, as equiv- 
alence transformations on a selected narrowed set of monitor functions. An obvious form that 
satisfies this requirement is the arc-length monitor function p = yl + au x . Alternatively, one 
could use, e.g., the similar functions p = yl + a ^x an d P = \/l + au x + (3h x or the curvature- 
related monitor functions p = \/l + au xx , p = \/l + ol\t? xx and p = \/l + au xx + (3h xx , where 
a and j3 are positive constants. 

Remark 2. The method of constructing an invariant discretization of a differential equation 
in combination with a numerical grid generator was discussed, e.g., in [12, 13, 16]. In contrast 
to the method employed above, in [21] the space of stencil variables was also prolonged to the 
monitor function, which is not necessarily chosen in an invariant way. In order to arrive at a 
completely invariant model, we however regard it important that the equidistribution principle 
is discretized in an invariant fashion too; see also the discussion in Section 6. Moreover, as the 
monitor function involves independent variables, unknown functions and their derivatives, it is 
possible to express its discretization using the same basis difference invariants of stencil variables 
that is needed for the physical differential equation discretization. In other words, no explicit 
prolongation to the monitor function is necessary within the framework of our approach. 

In Fig. 3 we show the integration of the one-dimensional shallow-water equations using the 
scheme (15), (17) with arc-length monitor function discretization (18) utilizing the same initial 
conditions as those chosen for the integration shown in Fig. 2. It is clearly visible that the mesh 
points remain almost fixed as long as the shock is not developed. Once the shock is traveling 
through the domain, the mesh points are able to sufficiently adapt to yield increased resolution 
in the region near the shock (as additionally shown in Fig. 3c). Again the scheme approximately 
conserves mass and momentum but dissipates energy. The relative errors in the momentum and 
energy conservation are approximately the same as in the case of the Lagrangian schemes in the 
previous subsection. 

It is worth pointing out that the time step of the integration shown in Fig. 3 is relatively 
small. The reason for this is that by using the scheme (15) and (17)-(18) we decouple the 
solution of the physical differential equation and the equation controlling the location of grid 
points. If time steps were not small, a severe time lag in the mesh movement would occur 
and the resulting mesh would not satisfy the equidistribution principle closely enough to give a 
satisfactory adaptivity. The above problem was extensively addressed in [29]. It can be overcome 
via the iterative solution of the physical and mesh equations a number of times, which leads to 
a reduction of the time lag in the mesh movement. Such a strategy could be readily adopted 
with the scheme (15) and (17)-(18) because a repeated iterative integration does not break the 
invariance of this scheme. 

In order to facilitate the comparison of distinct types of invariant numerical schemes, we also 
keep the time step small in the integration of the Lagrangian scheme for the one-dimensional 
shallow- water equations, which is presented in Fig. 2. 

5 Invariant numerical models for the two-dimensional 
shallow-water equations 

5.1 Selection of symmetries using boundary conditions 

The domains most often considered in geophysical fluid dynamics for the numerical integration 
of the two-dimensional shallow-water equations on a plane are either a channel with periodic 
boundary conditions in the East-West direction and rigid boundaries in the North-South di- 
rection or a domain with double periodic boundary conditions. As the second configuration 
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Figure 3. Numerical integration of the one-dimensional shallow-water equations (2) using the scheme (15) with 
r = 0.001 and N = 51 grid points on the domain [0, 2n] over the time interval [0,3]. The initial conditions are 
u = ylsiniE and h = ho + Asin(x + (po), with A — 0.4, </?o = 7r/6 and ho = 10. The trapezoidal rule is used for 
time integration and the arc-length monitor function is chosen for grid adaptation, setting a = 0.8. (a) Evolution 
of the discretization grid, (b) Numerical solution for h at t — 3. (c) Magnitude of the derivative u x of the solution 
for the scheme (15). Light shades refer to high values of \u x \. 



is more challenging from the point of view of invariant numerical schemes, we will employ it 
subsequently. 

Lie symmetry operators of the two-dimensional shallow-water equations (1) with periodic 
boundary conditions in both the East-West and North-South directions form the five-dimen- 
sional subalgebra S2 of the maximal Lie invariance algebra 02 of the equations (1) without 
additional constraints. A basis of S2 is given by 

d t , d x , d y , td x + d u , td y + d v . (19) 

As in Section 4 we could additionally include the scaling symmetries of the equations (1) in the 
subalgebra 52, referring to them as equivalence transformations of the class of doubly periodic 
boundary value problems. The reason why we did not include these scalings above is that 
all the discretizations for the shallow-water equations we use subsequently do not change the 
scaling properties of that system. This means that these discretizations already satisfy the 
required scaling properties by construction. On the other hand, the additional presence of scaling 
operators would slightly complicate the expressions for the difference invariants computed below, 
without giving any significant new information (the additional coefficients arising will factor out 
anyway for the resulting schemes). Only in the course of setting up the invariant grid generator 
will it be necessary to explicitly take into account the specific scaling symmetries, which will 
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consistently be done in Section 5.3. Note that both symmetries (19) and scaling symmetries of 
the two-dimensional shallow- water equations (1) generate equivalence transformations of the set 
of relevant initial conditions. 

As can be envisioned from the consideration of the numerical models of the one-dimensional 
shallow-water equations discussed in Section 4, discretization schemes for the two-dimensional 
shallow-water equations will be invariant under Galilean symmetries only if they are based on 
adaptive grids. Because for the channel model Galilean transformations are only admitted in the 
rr-direction, it suffices for a grid to be adaptive in the x-direction. This in particular means that 
we can use a uniform spacing in the y-direction and have a spatial grid with changing resolution 
only along the channel. On the other hand, the shallow-water equations with double periodic 
boundary conditions require the treatment of adaptive grids in both the x- and y-directions. 
An initial orthogonal spatial grid is driven to a nonorthogonal grid, which makes the direct 
evaluation of finite difference derivatives much more elaborate. This problem is treated upon 
using a finite volume formulation of that scheme. 

For simplicity, all the schemes are developed on an Arakawa A-grid subsequently; i.e., the 
variables u, v, h are defined in the same respective points. See, e.g., [42] for a discussion of 
different types of staggered grids for the shallow-water equations. Other types of grid staggering 
can be used in a similar way as that shown for the A-grid in what follows. 



5.2 Invariant numerical schemes with double periodic boundary conditions: 
Lagrangian scheme 

The main difficulty with adaptive grids in both the x- and y-directions is that it can become 
cumbersome to directly evaluate the gradients of the dependent variables on such curvilinear 
grids by finite differences. As discussed in Section 4, a prominent strategy to overcome this 
difficulty is to introduce a mapping from the computational (logical) coordinates (£, r]) to the 
physical coordinates (x,y). This will be done in Section 5.3. 

For the sake of demonstration we take another, more direct approach here, namely using 
the finite volume formulation of the divergence operator; see, e.g., [44]. Using the theorem of 
GauB-Ostrogradsky, we can approximate the divergence V • f of a vector-function f over a single 
grid cell with area A and edge lengths li as 

i 4 

V-f «-^(f i -n i )4. 

i=l 

In the above formula, rij denotes the outward-directed unit vector at the single cell edges. 
As it is possible to cast the shallow-water equations (1) into conserved (momentum) form, 



Eq h = h t + (hu) x + {hv) y = 0, 

Eq u = (hu) t + (hu 2 + hi 2 ^j + (huv) y = 0, 

Eq" = (hv) t + (huv) x + (hv 2 + -h 2 ) = 0, 



(20) 



2 



i) 



the above approximation of the divergence operator is sufficient to discretize the two-dimensional 
shallow- water equations using the finite volume form. 1 On the other hand, a finite volume 
discretization is readily applicable on adaptive grids, as it is not necessary to approximate 
derivatives by finite differences in such a formulation. 

In order to discretize (20) in an invariant way using the Dorodnitsyn method, we would 
need a set of difference invariants and would need to construct the discretization using these 



1 Equations including curl terms can be converted into finite volume representation using the Stokes theorem. 
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invariants as building blocks for the numerical scheme. The problem with this approach is the 
same as reported in the one-dimensional case. The Galilean transformation t = t, x = x + e±t, 
y = y + Bit, h = h, u = u + £\, v = v + £2 maps the system (20) to 

Eq h = Eq h , Eq u = Eq u + £iEq ft , Eq v = Eq v + £ 2 Eq h ; (21) 

i.e., it leads to a combination of the momentum equations with the continuity equation. Ex- 
pressing the momentum equations in terms of differential invariants thus again only works by 
combining these equations with the continuity equation. As in the one-dimensional case, it is 
therefore not natural to attempt to find an invariant approximation of the momentum form of 
the shallow- water equations using difference invariants. 

At the same time, the Lagrangian form of the shallow- water equations (1), which is 

dx dy dh , ( du dv\ du dh dv dh , . 

J="' di+Ka^r 0, * + * =0 - * + % =0 ' (22> 

can be approximated using the finite volume method as well. Expressing an invariant discretiza- 
tion of (22) in terms of difference invariants is considerably easier than doing the same for an 
invariant discretization of (20). As in the one-dimensional case, the drawback of using (22) as 
a starting point is that the resulting scheme does not approximate a conserved form and thus 
preserves neither mass and momenta nor energy. 

The stencil of the discretizations we aim to use is given in Fig. 4. All the dependent variables 
are defined in the centroids of the respective polygons. The fluxes through the edges will 
govern the evolution of these centroid values. In order to facilitate the computation of the 
fluxes it is necessary to determine the values of w = (u,v,h) in the cell corners, which is 
done by interpolation. While in principle any type of interpolation can be used, we employ 
natural neighbors interpolation for this purpose, i.e., the values at the cell corners are Wj = 
Yl K =i P K ^ w o' J 1 wnere 3 = 1, - - - , 4 and Wq' 3 are the values of w in the centers of those cells 
having in common the corner denoted by j. The interpolation weights p K ' 3 are determined 
in the following way. The Voronoi tessellation generated by the cell centers is constructed. 
Then a new tessellation is computed in which the point (xj,yj) is introduced as an additional 
generator. Denote by Ap* the area of the Voronoi cell of the original tessellation associated with 
the center point Pfi = (%o,yo) and by Ap j the area of the new cell associated with the corner 
point Pj = (xj,yj) introduced for the second tessellation. Then the weights p K ' 3 are computed 
as p K, i = {Ap 3 n Ap^)/Ap r Once the values Wj are obtained, they can be regarded as proper 
stencil variables. 

The prolongations of the symmetry operators (19) on the variables of the stencil shown in 
Fig. 4 read 



i=0 i=0 

4 4 
J2(td Xi + (t + r)d^ + d Ur + d &i ), J2(td yi + (t + r)% + d Vi + 8^). 



i=0 i=0 

These prolongations are well agreed with the above interpolation procedure. The difference 
invariants of the set (23) are given by 

t, hi, hi, Xi Xj f y% yj, Xi Xj tu^, yi yj tv^, 

Hi Uj , Ui Uj , Vi Vj , Vi Vj , 

where the indices i, j and k take the values 0, ... ,4. Note that of course not all of the above 
difference invariants are independent if i, j and k separately run through all possible values. 
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Figure 4. Stencil for the invariant Lagrangian schemes for the two-dimensional shallow-water equations with 
double periodic boundary conditions. The dependent variables w = (u,v,h) axe defined in the center (xo,J/o) of 
the cells of the respective polygons. The fluxes are computed using the values at the corners (xk,yk), k — 1 ... 4, 
of the cells, which are obtained by interpolation from the values at the polygon centroids. Variables with a hat 
are those at the subsequent time step. 



A simple explicit invariant scheme (Euler forward scheme) that can be constructed using 
these invariants is 

uq = 0, v = 0, 



X 




X 




T 




h 




h 




T 




Uo 




u 




T 








vo 


T 



+ 2^ ^2i( U i + u i+l)(Vi+l ~ Vi) - ( V i + V i+ i){x i+ i - Xi)\ = 0, 
i=l 

1 4 (24) 

i=l 
1 4 

i=l 

where A = \ Y,i=i( x iVi+l ~ x i+iVi) is the area of the polygon spanned by (xi, yi), . . . , (x 4 , y 4 ) 
and (X5, 2/5, u§, V5, h§) = (xi, y±, U\, V\, hi) by definition. As in the one-dimensional case, in the 
continuous limit this scheme converges to the Lagrangian representation of the two-dimensional 
shallow- water equations (22). 

In a similar manner, we can formulate the implicit scheme (trapezoidal rule) 

t;{ u o + uo) = 0, 7;i v o + vo) = 0, 

T l T I 

~ ~ + 77 ^2^ Ui + U i+l)(yi+l - yi) - ( V i + V i+ l)(x i+1 - Xi)] + 

i=l 

h 4 

TjX^C^ + "t+i)(z/i+l - &) - («i + v i+ i){x i+ i - Xi)) = 0, ^ 
iio - Uq , 1 



i=l 
4 



~~ + 77 + ~ Vi) + TT y](^i + hi+i)(yi+l - Vi) = 0, 

r 4/L ^— ' 44 



AA . 

i=i i=i 
1 



— — — - 77 y^Xhi + h i+ i)(x i+ i - Xi) - —j ^{hi + h i+1 )(x i+1 - £i) = 0. 



r f-f 4,4 f 
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Figure 5. Numerical integration of the two-dimensional shallow-water equations (1) using the scheme (25) with 
r = 0.001 and N x x N y = 71 x 71 grid points on the square [0, 2ir] x [0, 2ir] over the time interval [0, 2]. The initial 
conditions are u = A sin(x + tfio) sin y, v = A shirr sin y and h = ho + A cos (a; + ipo) cosy, with A = 0.4, tpo = 7r/6 
and ho = 10. Left: Numerical solution for h at t = 2. Right: Spatial discretization grid at t — 2. 

Fig. 5 shows the result of a numerical integration with the scheme (25) supplemented with 
periodic boundary conditions and specific initial conditions. The numerical solution of the 
water height h at t = 2 is shown in the left panel. The right panel depicts the associated 
discretization grid at t = 2. As in the case of the one-dimensional Lagrangian scheme, a 
strong distortion of the grid cells is visible, which is not directly related to pronounced features 
in the numerical solution, but rather a consequence of the Lagrangian grid movement. As 
both the discretizations (24) and (25) do not approximate the conserved form of the shallow- 
water equations (20), they conserve neither the mass A4 and the momenta V x and V y nor the 
energy H. 

It should be stressed that it is possible to formulate an invariant finite volume scheme for 
the conserved form of the shallow- water equations (20), in a similar way as was shown in the 
previous section for the one-dimensional case. As said above, the problem of doing this system- 
atically within the Dorodnitsyn approach is that it can be hard to find a proper combination 
of elementary difference invariants that allows one to approximate the momentum form of the 
shallow- water equations. This is why we will show an alternative way of constructing invari- 
ant numerical schemes for the two-dimensional shallow-water equations in the following section, 
which will avoid the technical complications that can arise when using the difference invariants 
method for the construction of symmetry-preserving numerical schemes. 

5.3 Invariant numerical schemes with double periodic boundary conditions: 
Eulerian scheme 

Though the finite volume discretization developed in Section 5.2 is suitable from the point of 
view of invar iance preservation, it is not ideal from the viewpoint of numerical analysis. In 
general, Lagrangian schemes are not in widespread use as they can easily lead to tangling 
meshes or rapidly changing grids through the spatial domain. For the same reason, numerical 
schemes in hydrodynamics are usually formulated in terms of Eulerian variables (or using some 
combination of Eulerian and Lagrangian schemes). An invariant scheme on an adaptive grid 
can be formulated by combining the idea of having an invariant grid generator proposed in 
Section 4.3 with the discretization in computational coordinates. More specifically, we consider 
the momentum form (20) of the two-dimensional shallow-water equations and rewrite it in the 
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computational coordinates 9 = t, £ = £(t,x,y), rj = r)(t, x, y): 



W9 {JFt) + lk {JitFt + JixFX + JiyFV) + lk {JlltFt + JVxFX + Jl]yFV) = °' (26) 
where F l = (h, hu, hv), F x = (hu, hu 2 + \h 2 , huv), F y = (hv, huv, hv 2 + \h?) and 

J = x £Vv ~ x vV^ & = ~ixxe - (.yVe, Vt = -Vx%6 - VyVe, 
a —VlL c - J^L - - 

The invariance of the system (26) with respect to shifts of the former independent variables t, 
x and y is obvious since the left-hand sides of equations of the system (which we denote by 
Eq fc , Eq" and Eq v , respectively) do not explicitly involve these variables. Therefore, any finite 
difference approximation of (26) is invariant with respect to the above shifts extended to the 
corresponding stencil; cf. (23). Note that any involved transformation is trivially extended to 
the computational coordinates £ and 77; i.e., they are not transformed. The scale symmetry 
transformations of the shallow- water equations (1) are automatically preserved in the course of 
a proper finite difference approximation of (26); see the related discussion in Section 5.1. In 
order to make clear the invariance with respect to Galilean boosts, we recombine terms in (26), 
substituting J£ t = -J£x%e ~ J^yVe and Jn t = -Jr] x x g - Jrj y y e : 

^{JVx(F x - xqF 1 ) + J Vy (F y - yeF 1 )) = 0. 

The Galilean boost t = t, x = x + £\t, y = y + e 2 t, h = h, u = u + £\, v = v + £2 maps the 
system (26) to the system with 

Eq /l = Eq ?t , Eq n = Eq n + eiEq h , Ekf = Eq v + e 2 Eq h . 

Note that this is the same transformation law in computational variables as it is in the physical 
space (21). The main idea for finding invariant numerical schemes of (26) is to construct the 
discretization in such a manner that the discrete counterpart of the system (26) is transformed 
similarly by the extension of the Galilean boost to the stencil points. In order to preserve Galilean 
boosts as symmetries in the course of discretization, it suffices 

• to use the same discretization schemes for all the equations of the system (26), just running 
through the number of corresponding components of F l , F x and F y ; 

• to evaluate all the copies of J£ x (resp. the components F and F x ) related to the block 
J£,x(F x — xqF 1 ) in the same grid point and in the same way; the same rule should be 
applied for the other similar blocks: J£y(F y — yoF 1 ), Ji] x (F x — xgF 1 ) and Jr] y (F y — yeF 1 ). 

For example, consider the trapezoidal scheme for the system (26): 

J Jk]jjk ~ J jkFjk + Ujk + Ujk + Vjk + Vjk = Q ^ 

where r is the step in 9 = t; 

Ujk = ■^^[(J&.)j+l/2,k( F jk + F j+l,k) ~ (J£t)j-l/2,k( F jk + F j-l,k) + 
(J£x)j+l/2,k{ F jk + F j+l,k) - (JCx)j-l/2,k( F jk + F j-l,k) + 
(Jd y ) j+ l/2,k( F Jk + F j + i,k) - (JQj-l/2,k( F Jk + F U,k)l 
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V jk = ^-[(JVt)j,k+i/2(F* k + Fl k+1 ) - (Jvt) j>k -i/2(Fjk + F j,k-i) + 
(JVx)j,k+i/2(F* k + F? )h+1 ) - {Jnx)j,k-i/2{Ff h + Ff >k _ x ) + 
{Jny)j,k+i/2{F] k + F£ k+1 ) - (J^) iffc _ 1/a (f* + FJ^)}; 

the values J, J£ x = y n , J£ y = —x n , Jr\ x = —y%, Jr\ y = x%, J£t and Jrj t are discretized in the 
following way: 

Jjk = K x j+l,k - x j~l,k)(Vj,k+l - Vj,k-l) - i x j,k+l ~ x j,k-l)(Vj+l,k - yj-l,k)], 

(J£x)j±i/2,k = -£^-(yj,k+i ~~ Vhk-l + Vj±i,k+i - Uj±i,k-i), 
{J£y)j±l/2,k = _ 4A^^'' fc+1 ~~ x h k ~ l + x j±l,k+l _ x j±l,k-l), 
(JVx)j,k±i/2 = —^^(Uj+l,k - Vj-l,k + Vj+l,k±i - Uj-l,k±l), 
{JVy)j,k±i/2 = 4^( x i+ 1 . fc ~~ + x j+i,k±i — x j-x,k±i), 

t it \ _ / it \ x jk + A J±i.fc ( jt \ + to.* 

W«Jj±i/2,fc — -W?xjj±i/2,fc 2 \ Jl iv)j±X/2,k 2 ' 

( T„ \ ( Tr, \ ± 3 k + ± i> k±1 ( TV, ^ y i k + y i' k±1 
[ Jr 1t)j,k±i/2 - -{ J Vx)j,k±i/2 2 K J Vy)j,k±l/2 ^ ' 

where Xj k = (£j k — Xj k )/r and yj k = (yj k — yj k )/r are by definition the mesh velocities in the x- 
and y-directions, respectively; and by hat we mark the corresponding values at the time 9 + r. 
In particular, we take 

rft\ -i ft \ x i k + x i± l > k i ft \ y i k + y i± l > k 

\ J ^t)j±i/2,k — \ Jt ix)j±\/2,k 2 \ J Sy)j±i/2,k 2 ' 

( Tm) - (fn) ±jk + ±j ' k±1 ( fn \ Vjk + Vj ' k±1 

\ jr lt)j,k±i/2 — \ jr lx)j tk ±i/2 2 \ jr ly)j,k±i/2 2 

As the system of difference equations (27) satisfies the above conditions, it is invariant with 
respect to properly extended Galilean boosts. 

Remark 3. The usage of computational coordinates also underlines the subtle change of the 
meaning of the time derivatives in a number of papers devoted to the construction of invariant 
numerical schemes, such as in [25, 50]. While in the standard (Eulerian) discretization, the 
continuous limit of the form (u — u)/t yields the partial derivative ut, in the framework of 
invariant schemes these terms are often to be interpreted as Lagrangian time derivatives u 
(see also Section 4.2). This immediate transition from an Eulerian (partial) time derivative 
to a Lagrangian (total) time derivative is a necessary consequence of the intermediate step of 
discretizing an equation in computational coordinates and assuming that the grid evolution is 
described by the equations x = u and y = v; see also the discussion in [43]. 

Remark 4. It should be noted that computational coordinates have a clear physical meaning 
in the present context. As they do not change during the evolution of the grid, they can be 
interpreted as the Lagrangian variables (fluid labels) of fluid mechanics provided we again assume 
a Lagrangian grid evolution. A prominent way to choose these fluid labels is by setting them to 
equal the Cartesian coordinates at the onset of evolution. By definition, this is the same role 
that computational coordinates play in the numerics of moving meshes. Stated in another way, 
the requirement of maintaining invariance of the discretization scheme and discretization stencil 
of the shallow-water equations under the Galilean group naturally boils down to discretizing 
these equations in computational coordinates. 
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It then remains to specify the grid velocities x and y in order to complete the scheme given 
in (27). This can be done in a similar manner as in Section 4.3 using the idea of equidistributing 
meshes (though, strictly speaking, equidistribution in higher dimensions is not sufficient to 
uniquely determine an adaptive grid; see, e.g., the discussion in [29]). Thus, the grid will be 
determined from the system of elliptic equations 

V € • (GV c x) =0, • (GV^y) = 0. (28) 

Here denotes the gradient in the space of computational coordinates (£, rj) and G = wl is 
the matrix-valued monitor function, where I is the two-by-two unit matrix and w = w(x, y) is a 
weight function which depends on the (numerical) solution of the shallow- water equations [18]. 
An invariant discretization of (28) reads 

1 / %ij z ij z i—l,j 



1 l z i,j+l z ij z ij z i,j-l 1 

Wi_n>i o : Wi o ; — U 



Arj ^ "ij-i/2- Ar] 

where z = x and z = y for the first and second equations, respectively, and 



w i+l/2,j = 7. > w i-l/2,3 



Wi + ij+ Wij _ _ Wij + 

2 

Wij+l + Wij Wij + Wij-l 

u>ij+i/2 = 2 ' Wi >i-V 2 = 2 

provided that w is approximated by a difference invariant of the algebra 52, which is spanned by 
the vector fields (19). Once again, straightforward choices for w that can be discretized using 
difference invariants are 



w 1 = J I + a{ul + Uy\ + vl+ v%), w 2 = J 1 + a(h xx + h yy ) 2 , 



but of course other forms for w are possible as well. Similarly to the one-dimensional case (cf. 
Remark 1), the general form of w that is an invariant of the algebra S2 is given by an arbitrary 
smooth function of derivatives of u, v and h with respect to x and y including h itself but not u 
and v. At the same time, we should also take into account other desired properties of w. Both 
the functions w 1 and w 2 are invariant with respect to shifts and Galilean boosts generated by 
vector fields (19), the scalings generated by the vector field xd x + yd y + ud u + vd v + 2hdh and 
even rotations. All the scaling symmetries of the shallow- water equations are at least equivalence 
transformations for the sets of functions of such forms, where the parameter a is varied. 

It should also be stressed that the grid generator based on system (28) is a rather simple one. 
More advanced formulations of grid generators exist, e.g., by using a general positive definite 
and symmetric matrix G. Alternatively, the grids at a certain time level could be computed 
using moving mesh partial differential equations [15, 29], provided it is possible to discretize 
such equations in an invariant way. 

A different methodology is to use so-called velocity-based moving mesh strategies. Unlike 
the location-based methods, which were exclusively used in the present paper, in the velocity- 
based methods the location of the grid points is not determined directly but rather equations 
for the mesh velocity are formulated. Velocity based strategies, such as the method involving 
the geometric conservation law [17, 29], provide alternative ways of formulating grid equations 
that give a basis for realizing invariant moving mesh equations. See also [48], where the term 
geometric conservation law was introduced. 

In Fig. 6 we repeat the numerical integration of the two-dimensional shallow-water equations 
with the setting of Fig. 5 but now using the scheme (27) in combination with a grid generator 
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Figure 6. Numerical integration of the two-dimensional shallow-water equations (1) using the scheme (27) with 
r = 0.001 and N x x N y = 71 x 71 grid points on the square [0, 2ir] x [0, 2ir] over the time interval [0, 2]. The initial 
conditions are u = A sm(x + (fio) siny, v = A sin x sin y and h — ho + A cos (a; + fo) cos y, with A = 0.4, ipo = 7r/6 
and ho = 10. As a weight function, w 2 = ^/l + a(h xx + h yy ) 2 is chosen with a = 0.4. (a) Numerical solution for 
h at £ = 2. (6) Spatial discretization grid at t = 2. (c) The weight function at t = 2. 



based on w 2 . Similarly to the case of the one-dimensional shallow- water equations, it can be 
seen from Fig. 6 that the usage of a grid generator leads to grids that are not as distorted as 
those obtained from a purely Lagrangian scheme. Moreover, the regions of grid concentration 
are now directly linked to the physical behavior of the numerical solution for the dependent 
variable h. The scheme (27) is mass and momentum conserving but, as all the other schemes 
presented in the paper, dissipates the energy. The conserved quantities are evaluated at time t 
as M = A^Ar]Y / j i k h jkJjk, = A^Ar]J2j,k h jkUjkJjk, Vy = A^Ar]Y / j i k h jkVjkJjk and % = 

6 Conclusion 

The present paper is devoted to the construction of several invariant numerical schemes mod- 
eling shallow-water dynamics. In particular, we aim to describe a possible bridge between the 
formalism of constructing invariant numerical schemes and the existing theory on adaptive mov- 
ing meshes. Such a bridge was already indicated in the literature. Indeed, there already exist a 
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number of investigations devoted to the importance of scale invariance in the theory of moving 
mesh equations. Thus, in [14, 15, 16] (see also [29, p. Ill] and references therein) a moving mesh 
partial differential equation was constructed that preserves the scale invariance of the physical 
differential equation to be discretized. The extension of this idea to setting up a grid generator 
that is invariant with respect to (a suitable subgroup of) the maximal Lie invariance group of 
a system of differential equations is therefore straightforward and was conceptually indicated 
in the aforementioned sources. The idea of introducing computational coordinates into invari- 
ant numerical schemes has also been successfully demonstrated for one-dimensional nonlinear 
Schrodinger equations [11]. 

We require our discretization schemes to be invariant with respect to the subgroup of the 
maximal Lie symmetry group of the (resp. one- or two-dimensional) shallow- water equations ad- 
mitted when imposing periodic boundary conditions. From the physical point of view it is natural 
to assume that appropriate symmetries of the system of differential equations under consider- 
ation act as equivalence transformations on a joint class of physically relevant boundary value 
problems. Imposing periodic boundary conditions for varying intervals in the one-dimensional 
case (resp. for rectangular domains of varying sizes whose sides are parallel to coordinate axes 
in the two-dimensional case) while both the initial time and initial conditions also vary, we 
select the subgroup generated by the time and space translations, the Galilean boosts and the 
scalings symmetries of the shallow-water equations. Other subgroups might be chosen as well, 
but for wide or even infinite-dimensional maximal Lie invariance (pseudo)groups admitted by 
the prominent models in hydrodynamics it can be quite intricate to justify the choice for such 
subgroups from the physical point of view. 

In general, the inclusion of well-proven principles in the study of invariant numerical schemes 
is required. The invariant schemes for numerous evolution equations constructed so far were 
mostly purely Lagrangian schemes. However, these schemes are not in prevalent use in practice as 
they usually lead to complicated mesh geometries that might eventually (at least locally) degrade 
the quality of the numerical solution. This can be seen directly by comparing Figures 5 and 6, 
where the stronger distortion of the grid lines in the Lagrangian scheme is already manifest after 
a relatively short period of integration. Therefore, the formulation of invariant grid generators 
coupled with suitable invariant discrete counterparts of physical systems of differential equations 
is a practical way for symmetry preserving numerical integration of these systems. 

A further novel feature of the present paper is the construction of invariant numerical schemes 
for higher-dimensional systems of partial differential equations. Higher-dimensional schemes 
are especially challenging if it is not possible to use fixed orthogonal grids. In the course 
of constructing such schemes for the two-dimensional shallow-water equations we have shown 
that invariant discretizations are not only restricted to finite difference schemes. It is possible 
and straightforward to also formulate finite volume discretizations that preserve symmetries of 
systems of differential equations. In a similar manner, other discretization techniques, such as 
the finite element method, could be employed as well. 

This problem can also be tackled by transforming the system under consideration into com- 
putational coordinates. We have shown that the transition to these coordinates is a natural step 
in the course of the construction of Galilean invariant discrete schemes. The key to the construc- 
tion is then not to simply combine difference invariants as proposed in the original method by 
Dorodnitsyn but to study the transformation laws of the equations in computational coordinates 
for the respective symmetries. These laws are trivial for all the admitted symmetries except for 
the Galilean boosts. For Galilean invariance it is found that the new momentum equations are 
given as the combination of the old momentum and continuity equations. An invariant dis- 
cretization is therefore achieved by finding proper discrete counterparts of these transformation 
laws rather than combining difference invariants. 

Because the main objective of this paper is to demonstrate different strategies for finding 
invariant discretization schemes exemplified with the shallow- water equations, the discretization 
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schemes considered are kept as simple as possible. This concerns both the design of the schemes 
themselves (e.g., using only two-level time integration methods and unstaggered grids) and the 
solution of the resulting algebraic equations, which is done in the most direct and straightforward 
manner. More advanced integration and algebraic solution techniques can be readily adopted 
but their discussion is restrained to keep the focus of the paper on the conceptual aspect of 
introducing the Lie symmetry approach in the framework of numerical modeling as far as pos- 
sible. For example, the extension to more advanced time integration methods such as arbitrary 
Runge-Kutta and general time-splitting schemes can be done by extension of the discretization 
stencils via inclusion of further time layers. Similarly, the usage of staggered grids can be facili- 
tated by adding further points to the stencil on the same time layer. The procedure of invariant 
discretization involving wider stencils then follows precisely the same techniques as outlined and 
used in the present paper. Within the approach based on the construction of difference invari- 
ants, both ways of extending the stencils will lead to a larger number of invariants and thus to an 
increased number of possibilities for combining them to form a particular discretization scheme. 

It was mentioned in the introduction that a system of differential equations might possess 
various qualitative properties that one should aim to preserve in the course of setting up a nu- 
merical model. Besides symmetries, it is of outstanding importance to monitor the behavior of 
conserved quantities possessed by the system under consideration. This is a problem of cen- 
tral importance in long-term integrations of such systems as a systematic failure in capturing 
conservation laws may lead to unrealistic numerical results (e.g., loss of mass or wrong turbu- 
lence spectra). Proper discretizations of the momentum form of the shallow- water equations 
conserve the mass and momenta exactly or to high order, but none of them is actually energy 
conserving. This should not come as a complete surprise as setting up energy-conserving schemes 
for the shallow- water equations is a quite nontrivial problem; see, e.g., the schemes proposed 
in [2, 45, 46]. The inclusion of additional conserved quantities in the construction of invariant 
discretization schemes will therefore be one of our future research topics. 

From a more general point of view, the requirement of preserving symmetries in a discretiza- 
tion scheme might lead to a geometric justification for using adaptive meshes. Though there 
are several classes of physical problems (such as blow-ups) for which adaptive meshes are well 
suited, the usage of such meshes is not undisputed in the numerical analysis and geophysical 
fluid dynamics communities. The drawbacks of moving meshes, such as an additional level of 
complexity of the schemes and the computational overhead resulting from computing and stor- 
ing the mesh points at each time level, must be well opposed to their potential benefits on a 
case-by-case basis. The result that the numerical preservation of important structural properties 
like symmetries automatically requires using moving meshes can thus be seen as a geometric 
argument for allowing adaptive discretization grids for certain classes of physical differential 
equations. Moreover, the usage of a grid redistribution equation (or r-adaptivity) as advocated 
in the present paper is also most suitable because it can be efficiently implemented within the 
framework of parallel computing. 
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